import xarray as xr
import xesmf as xe
import hvplot.xarray
import pandas as pd
import os
import glob
import geopandas as gpd
from distributed import Client
import matplotlib.pyplot as plt
os.chdir('/Users/caramelo/Documents/GitHub/cmc/')
grib_list_total=sorted(glob.glob('test_reg*.grib2'))
#grib_list=grib_list_total[0:7]
#grib_list
grib_list_total
['test_reg202103260.grib2', 'test_reg202103261.grib2', 'test_reg2021032610.grib2', 'test_reg2021032611.grib2', 'test_reg2021032612.grib2', 'test_reg2021032613.grib2', 'test_reg2021032614.grib2', 'test_reg2021032615.grib2', 'test_reg2021032616.grib2', 'test_reg2021032617.grib2', 'test_reg2021032618.grib2', 'test_reg2021032619.grib2', 'test_reg202103262.grib2', 'test_reg2021032620.grib2', 'test_reg2021032621.grib2', 'test_reg2021032622.grib2', 'test_reg2021032623.grib2', 'test_reg2021032624.grib2', 'test_reg2021032625.grib2', 'test_reg2021032626.grib2', 'test_reg2021032627.grib2', 'test_reg2021032628.grib2', 'test_reg2021032629.grib2', 'test_reg202103263.grib2', 'test_reg2021032630.grib2', 'test_reg2021032631.grib2', 'test_reg2021032632.grib2', 'test_reg2021032633.grib2', 'test_reg2021032634.grib2', 'test_reg2021032635.grib2', 'test_reg2021032636.grib2', 'test_reg2021032637.grib2', 'test_reg2021032638.grib2', 'test_reg2021032639.grib2', 'test_reg202103264.grib2', 'test_reg2021032640.grib2', 'test_reg2021032641.grib2', 'test_reg2021032642.grib2', 'test_reg2021032643.grib2', 'test_reg2021032644.grib2', 'test_reg2021032645.grib2', 'test_reg2021032646.grib2', 'test_reg2021032647.grib2', 'test_reg2021032648.grib2', 'test_reg2021032649.grib2', 'test_reg202103265.grib2', 'test_reg2021032650.grib2', 'test_reg2021032651.grib2', 'test_reg2021032652.grib2', 'test_reg2021032653.grib2', 'test_reg2021032654.grib2', 'test_reg2021032655.grib2', 'test_reg2021032656.grib2', 'test_reg2021032657.grib2', 'test_reg2021032658.grib2', 'test_reg2021032659.grib2', 'test_reg202103266.grib2', 'test_reg2021032660.grib2', 'test_reg2021032661.grib2', 'test_reg2021032662.grib2', 'test_reg2021032663.grib2', 'test_reg2021032664.grib2', 'test_reg2021032665.grib2', 'test_reg2021032666.grib2', 'test_reg2021032667.grib2', 'test_reg2021032668.grib2', 'test_reg2021032669.grib2', 'test_reg202103267.grib2', 'test_reg2021032670.grib2', 'test_reg2021032671.grib2', 'test_reg2021032672.grib2', 'test_reg2021032673.grib2', 'test_reg2021032674.grib2', 'test_reg2021032675.grib2', 'test_reg2021032676.grib2', 'test_reg2021032677.grib2', 'test_reg2021032678.grib2', 'test_reg2021032679.grib2', 'test_reg202103268.grib2', 'test_reg2021032680.grib2', 'test_reg2021032681.grib2', 'test_reg2021032682.grib2', 'test_reg2021032683.grib2', 'test_reg2021032684.grib2', 'test_reg202103269.grib2']
ds=xr.open_mfdataset(grib_list_total,concat_dim='valid_time',engine='cfgrib',combine='nested')
ds
<xarray.Dataset>
Dimensions: (valid_time: 85, x: 935, y: 824)
Coordinates:
time datetime64[ns] 2021-03-26
step (valid_time) timedelta64[ns] 00:00:00 ... 0 days 09:00:00
heightAboveGround int64 10
latitude (y, x) float64 dask.array<chunksize=(824, 935), meta=np.ndarray>
longitude (y, x) float64 dask.array<chunksize=(824, 935), meta=np.ndarray>
* valid_time (valid_time) datetime64[ns] 2021-03-26 ... 2021-03-26T...
Dimensions without coordinates: x, y
Data variables:
gust (valid_time, y, x) float32 dask.array<chunksize=(1, 824, 935), meta=np.ndarray>
Attributes:
GRIB_edition: 2
GRIB_centre: cwao
GRIB_centreDescription: Canadian Meteorological Service - Montreal
GRIB_subCentre: 0
Conventions: CF-1.7
institution: Canadian Meteorological Service - Montreal
history: 2021-03-26T13:58:15 GRIB to CDM+CF via cfgrib-0....array('2021-03-26T00:00:00.000000000', dtype='datetime64[ns]')array([ 0, 3600000000000, 36000000000000, 39600000000000,
43200000000000, 46800000000000, 50400000000000, 54000000000000,
57600000000000, 61200000000000, 64800000000000, 68400000000000,
7200000000000, 72000000000000, 75600000000000, 79200000000000,
82800000000000, 86400000000000, 90000000000000, 93600000000000,
97200000000000, 100800000000000, 104400000000000, 10800000000000,
108000000000000, 111600000000000, 115200000000000, 118800000000000,
122400000000000, 126000000000000, 129600000000000, 133200000000000,
136800000000000, 140400000000000, 14400000000000, 144000000000000,
147600000000000, 151200000000000, 154800000000000, 158400000000000,
162000000000000, 165600000000000, 169200000000000, 172800000000000,
176400000000000, 18000000000000, 180000000000000, 183600000000000,
187200000000000, 190800000000000, 194400000000000, 198000000000000,
201600000000000, 205200000000000, 208800000000000, 212400000000000,
21600000000000, 216000000000000, 219600000000000, 223200000000000,
226800000000000, 230400000000000, 234000000000000, 237600000000000,
241200000000000, 244800000000000, 248400000000000, 25200000000000,
252000000000000, 255600000000000, 259200000000000, 262800000000000,
266400000000000, 270000000000000, 273600000000000, 277200000000000,
280800000000000, 284400000000000, 28800000000000, 288000000000000,
291600000000000, 295200000000000, 298800000000000, 302400000000000,
32400000000000], dtype='timedelta64[ns]')array(10)
|
|
array(['2021-03-26T00:00:00.000000000', '2021-03-26T01:00:00.000000000',
'2021-03-26T10:00:00.000000000', '2021-03-26T11:00:00.000000000',
'2021-03-26T12:00:00.000000000', '2021-03-26T13:00:00.000000000',
'2021-03-26T14:00:00.000000000', '2021-03-26T15:00:00.000000000',
'2021-03-26T16:00:00.000000000', '2021-03-26T17:00:00.000000000',
'2021-03-26T18:00:00.000000000', '2021-03-26T19:00:00.000000000',
'2021-03-26T02:00:00.000000000', '2021-03-26T20:00:00.000000000',
'2021-03-26T21:00:00.000000000', '2021-03-26T22:00:00.000000000',
'2021-03-26T23:00:00.000000000', '2021-03-27T00:00:00.000000000',
'2021-03-27T01:00:00.000000000', '2021-03-27T02:00:00.000000000',
'2021-03-27T03:00:00.000000000', '2021-03-27T04:00:00.000000000',
'2021-03-27T05:00:00.000000000', '2021-03-26T03:00:00.000000000',
'2021-03-27T06:00:00.000000000', '2021-03-27T07:00:00.000000000',
'2021-03-27T08:00:00.000000000', '2021-03-27T09:00:00.000000000',
'2021-03-27T10:00:00.000000000', '2021-03-27T11:00:00.000000000',
'2021-03-27T12:00:00.000000000', '2021-03-27T13:00:00.000000000',
'2021-03-27T14:00:00.000000000', '2021-03-27T15:00:00.000000000',
'2021-03-26T04:00:00.000000000', '2021-03-27T16:00:00.000000000',
'2021-03-27T17:00:00.000000000', '2021-03-27T18:00:00.000000000',
'2021-03-27T19:00:00.000000000', '2021-03-27T20:00:00.000000000',
'2021-03-27T21:00:00.000000000', '2021-03-27T22:00:00.000000000',
'2021-03-27T23:00:00.000000000', '2021-03-28T00:00:00.000000000',
'2021-03-28T01:00:00.000000000', '2021-03-26T05:00:00.000000000',
'2021-03-28T02:00:00.000000000', '2021-03-28T03:00:00.000000000',
'2021-03-28T04:00:00.000000000', '2021-03-28T05:00:00.000000000',
'2021-03-28T06:00:00.000000000', '2021-03-28T07:00:00.000000000',
'2021-03-28T08:00:00.000000000', '2021-03-28T09:00:00.000000000',
'2021-03-28T10:00:00.000000000', '2021-03-28T11:00:00.000000000',
'2021-03-26T06:00:00.000000000', '2021-03-28T12:00:00.000000000',
'2021-03-28T13:00:00.000000000', '2021-03-28T14:00:00.000000000',
'2021-03-28T15:00:00.000000000', '2021-03-28T16:00:00.000000000',
'2021-03-28T17:00:00.000000000', '2021-03-28T18:00:00.000000000',
'2021-03-28T19:00:00.000000000', '2021-03-28T20:00:00.000000000',
'2021-03-28T21:00:00.000000000', '2021-03-26T07:00:00.000000000',
'2021-03-28T22:00:00.000000000', '2021-03-28T23:00:00.000000000',
'2021-03-29T00:00:00.000000000', '2021-03-29T01:00:00.000000000',
'2021-03-29T02:00:00.000000000', '2021-03-29T03:00:00.000000000',
'2021-03-29T04:00:00.000000000', '2021-03-29T05:00:00.000000000',
'2021-03-29T06:00:00.000000000', '2021-03-29T07:00:00.000000000',
'2021-03-26T08:00:00.000000000', '2021-03-29T08:00:00.000000000',
'2021-03-29T09:00:00.000000000', '2021-03-29T10:00:00.000000000',
'2021-03-29T11:00:00.000000000', '2021-03-29T12:00:00.000000000',
'2021-03-26T09:00:00.000000000'], dtype='datetime64[ns]')
|
ds.gust.attrs
{'GRIB_paramId': 260065,
'GRIB_shortName': 'gust',
'GRIB_units': 'm s**-1',
'GRIB_name': 'Wind speed (gust)',
'GRIB_cfVarName': 'gust',
'GRIB_dataType': 'af',
'GRIB_missingValue': 9999,
'GRIB_numberOfPoints': 770440,
'GRIB_typeOfLevel': 'heightAboveGround',
'GRIB_NV': 0,
'GRIB_stepUnits': 1,
'GRIB_stepType': 'instant',
'GRIB_gridType': 'polar_stereographic',
'GRIB_gridDefinitionDescription': 'Polar stereographic ',
'long_name': 'Wind speed (gust)',
'units': 'm s**-1'}
ds.gust.isel(valid_time=1).plot()
<matplotlib.collections.QuadMesh at 0x167e8c9a0>
ds = ds.rename({'latitude': 'lat', 'longitude': 'lon'})
#ds
ds_out = xe.util.grid_global(0.1, 0.1)
#ds_out # contains lat/lon values of cell centers and boundaries.
regridder = xe.Regridder(ds, ds_out, 'bilinear')
regridder
dr_out = regridder(ds)
dr_out.where(dr_out>0).gust.plot()
dr_out['lat'] = dr_out.lat[:,0].drop('lon')
dr_out['lon'] = dr_out.lon[0,:].drop('lat')
dr_out = dr_out.sortby(['x','y'])
dr_out = dr_out.rename({'x':'longitude',
'y':'latitude',
'lon':'longitude',
'lat':'latitude'})
dr_out = dr_out.where(dr_out>0) # remplacer 0 par nan
/Users/caramelo/anaconda3/envs/test_xesmf/lib/python3.9/site-packages/xesmf/util.py:94: UserWarning: 360 cannot be divided by d_lon = 0.1, might not cover the globe uniformally warnings.warn( /Users/caramelo/anaconda3/envs/test_xesmf/lib/python3.9/site-packages/xesmf/util.py:100: UserWarning: 180 cannot be divided by d_lat = 0.1, might not cover the globe uniformally warnings.warn( /Users/caramelo/anaconda3/envs/test_xesmf/lib/python3.9/site-packages/xesmf/frontend.py:487: FutureWarning: ``output_sizes`` should be given in the ``dask_gufunc_kwargs`` parameter. It will be removed as direct parameter in a future version. ds_out = xr.apply_ufunc(
dr_out
<xarray.Dataset>
Dimensions: (latitude: 1800, longitude: 3600, valid_time: 85)
Coordinates:
time datetime64[ns] 2021-03-26
step (valid_time) timedelta64[ns] 00:00:00 ... 0 days 09:00:00
heightAboveGround int64 10
* valid_time (valid_time) datetime64[ns] 2021-03-26 ... 2021-03-26T...
* longitude (longitude) float64 -179.9 -179.9 -179.8 ... 179.8 179.9
* latitude (latitude) float64 -89.95 -89.85 -89.75 ... 89.85 89.95
Data variables:
gust (valid_time, latitude, longitude) float64 dask.array<chunksize=(1, 1800, 3600), meta=np.ndarray>
Attributes:
regrid_method: bilineararray('2021-03-26T00:00:00.000000000', dtype='datetime64[ns]')array([ 0, 3600000000000, 36000000000000, 39600000000000,
43200000000000, 46800000000000, 50400000000000, 54000000000000,
57600000000000, 61200000000000, 64800000000000, 68400000000000,
7200000000000, 72000000000000, 75600000000000, 79200000000000,
82800000000000, 86400000000000, 90000000000000, 93600000000000,
97200000000000, 100800000000000, 104400000000000, 10800000000000,
108000000000000, 111600000000000, 115200000000000, 118800000000000,
122400000000000, 126000000000000, 129600000000000, 133200000000000,
136800000000000, 140400000000000, 14400000000000, 144000000000000,
147600000000000, 151200000000000, 154800000000000, 158400000000000,
162000000000000, 165600000000000, 169200000000000, 172800000000000,
176400000000000, 18000000000000, 180000000000000, 183600000000000,
187200000000000, 190800000000000, 194400000000000, 198000000000000,
201600000000000, 205200000000000, 208800000000000, 212400000000000,
21600000000000, 216000000000000, 219600000000000, 223200000000000,
226800000000000, 230400000000000, 234000000000000, 237600000000000,
241200000000000, 244800000000000, 248400000000000, 25200000000000,
252000000000000, 255600000000000, 259200000000000, 262800000000000,
266400000000000, 270000000000000, 273600000000000, 277200000000000,
280800000000000, 284400000000000, 28800000000000, 288000000000000,
291600000000000, 295200000000000, 298800000000000, 302400000000000,
32400000000000], dtype='timedelta64[ns]')array(10)
array(['2021-03-26T00:00:00.000000000', '2021-03-26T01:00:00.000000000',
'2021-03-26T10:00:00.000000000', '2021-03-26T11:00:00.000000000',
'2021-03-26T12:00:00.000000000', '2021-03-26T13:00:00.000000000',
'2021-03-26T14:00:00.000000000', '2021-03-26T15:00:00.000000000',
'2021-03-26T16:00:00.000000000', '2021-03-26T17:00:00.000000000',
'2021-03-26T18:00:00.000000000', '2021-03-26T19:00:00.000000000',
'2021-03-26T02:00:00.000000000', '2021-03-26T20:00:00.000000000',
'2021-03-26T21:00:00.000000000', '2021-03-26T22:00:00.000000000',
'2021-03-26T23:00:00.000000000', '2021-03-27T00:00:00.000000000',
'2021-03-27T01:00:00.000000000', '2021-03-27T02:00:00.000000000',
'2021-03-27T03:00:00.000000000', '2021-03-27T04:00:00.000000000',
'2021-03-27T05:00:00.000000000', '2021-03-26T03:00:00.000000000',
'2021-03-27T06:00:00.000000000', '2021-03-27T07:00:00.000000000',
'2021-03-27T08:00:00.000000000', '2021-03-27T09:00:00.000000000',
'2021-03-27T10:00:00.000000000', '2021-03-27T11:00:00.000000000',
'2021-03-27T12:00:00.000000000', '2021-03-27T13:00:00.000000000',
'2021-03-27T14:00:00.000000000', '2021-03-27T15:00:00.000000000',
'2021-03-26T04:00:00.000000000', '2021-03-27T16:00:00.000000000',
'2021-03-27T17:00:00.000000000', '2021-03-27T18:00:00.000000000',
'2021-03-27T19:00:00.000000000', '2021-03-27T20:00:00.000000000',
'2021-03-27T21:00:00.000000000', '2021-03-27T22:00:00.000000000',
'2021-03-27T23:00:00.000000000', '2021-03-28T00:00:00.000000000',
'2021-03-28T01:00:00.000000000', '2021-03-26T05:00:00.000000000',
'2021-03-28T02:00:00.000000000', '2021-03-28T03:00:00.000000000',
'2021-03-28T04:00:00.000000000', '2021-03-28T05:00:00.000000000',
'2021-03-28T06:00:00.000000000', '2021-03-28T07:00:00.000000000',
'2021-03-28T08:00:00.000000000', '2021-03-28T09:00:00.000000000',
'2021-03-28T10:00:00.000000000', '2021-03-28T11:00:00.000000000',
'2021-03-26T06:00:00.000000000', '2021-03-28T12:00:00.000000000',
'2021-03-28T13:00:00.000000000', '2021-03-28T14:00:00.000000000',
'2021-03-28T15:00:00.000000000', '2021-03-28T16:00:00.000000000',
'2021-03-28T17:00:00.000000000', '2021-03-28T18:00:00.000000000',
'2021-03-28T19:00:00.000000000', '2021-03-28T20:00:00.000000000',
'2021-03-28T21:00:00.000000000', '2021-03-26T07:00:00.000000000',
'2021-03-28T22:00:00.000000000', '2021-03-28T23:00:00.000000000',
'2021-03-29T00:00:00.000000000', '2021-03-29T01:00:00.000000000',
'2021-03-29T02:00:00.000000000', '2021-03-29T03:00:00.000000000',
'2021-03-29T04:00:00.000000000', '2021-03-29T05:00:00.000000000',
'2021-03-29T06:00:00.000000000', '2021-03-29T07:00:00.000000000',
'2021-03-26T08:00:00.000000000', '2021-03-29T08:00:00.000000000',
'2021-03-29T09:00:00.000000000', '2021-03-29T10:00:00.000000000',
'2021-03-29T11:00:00.000000000', '2021-03-29T12:00:00.000000000',
'2021-03-26T09:00:00.000000000'], dtype='datetime64[ns]')array([-179.95, -179.85, -179.75, ..., 179.75, 179.85, 179.95])
array([-89.95, -89.85, -89.75, ..., 89.75, 89.85, 89.95])
|
os.chdir('/Users/caramelo/Documents/GitHub/cmc/')
coords_stations=pd.read_csv('Stations matrice scribe.csv')
coords_stations
| OACI | NOMSTN | LAT | LON | |
|---|---|---|---|---|
| 0 | CYUL | MONTREAL (DORVAL) | 45.4667 | -73.7500 |
| 1 | CYOW | OTTAWA | 45.3167 | -75.6667 |
| 2 | CYVO | VAL D'OR | 48.0500 | -77.7833 |
| 3 | CYSC | SHERBROOKE | 45.4333 | -71.6833 |
| 4 | CYQB | QUEBEC | 46.8000 | -71.4000 |
| 5 | CYRJ | ROBERVAL | 48.5167 | -72.2667 |
| 6 | CYGL | LA GRANDE RIVIERE | 53.6333 | -77.7000 |
| 7 | CYBC | BAIE COMEAU | 49.1333 | -68.2000 |
| 8 | CYZV | SEPT-ILES | 50.2167 | -66.2500 |
| 9 | CYRQ | TROIS-RIVIÈRES (A) | 46.3539 | -72.6792 |
| 10 | CYMX | MIRABEL | 45.6833 | -74.0333 |
| 11 | CWPK | PARENT | 47.9167 | -74.6167 |
| 12 | CYYY | MONT JOLI | 48.6167 | -68.2167 |
| 13 | CYWK | WABUSH LAKE | 52.9333 | -66.8667 |
| 14 | CYAH | LA GRANDE IV | 53.7500 | -73.6667 |
| 15 | CYYB | NORTH BAY | 46.3667 | -79.4167 |
| 16 | CYGP | GASPE | 48.7833 | -64.4833 |
#Montréal
gust1d = dr_out.gust.sel(latitude=coords_stations.iloc[0,2], longitude=coords_stations.iloc[0,3],method='nearest')
df = gust1d.reset_coords(drop=True).to_dataframe()
df.plot()
plt.title(coords_stations.iloc[0,1])
Text(0.5, 1.0, 'MONTREAL (DORVAL)')
#bon maintenant admettons pour la station de sherbrooke
gust1d = dr_out.gust.sel(latitude=coords_stations.iloc[3,2], longitude=coords_stations.iloc[3,3],method='nearest')
df = gust1d.reset_coords(drop=True).to_dataframe()
df.plot()
plt.title(coords_stations.iloc[3,1])
Text(0.5, 1.0, 'SHERBROOKE')
import holoviews as hv
from holoviews import opts
hv.extension('bokeh')
import hvplot.pandas # noqa
df.hvplot(title='Sherbrooke')
import hvplot.pandas # noqa
df.hvplot(title='Montréal')